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Abstract We study a class of shear-free, homogeneous but anisotropic cos¬ 
mological models with imperfect matter sources in the context of f{R) gravity. 
We show that the anisotropic stresses are related to the electric part of the 
Weyl tensor in such a way that they balance each other. We also show that 
within the class of orthogonal f{R) models, small perturbations of shear are 
damped, and that the electric part of the Weyl tensor and the anisotropic 
stress tensor decay with the expansion as well as the heat flux of the curvature 
fluid. Specializing in locally rotationally symmetric spacetimes in orthonormal 
frames, we examine the late-time behaviour of the de Sitter universe in f{R) 
gravity. For the Starobinsky model of f{R), we study the evolutionary behav¬ 
ior of the Universe by numerically integrating the Friedmann equation, where 
the initial conditions for the expansion, acceleration and jerk parameters are 
taken from observational data. 
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1 Introduction 

As a result of the current understanding that the Universe is in a state of accel¬ 
erated expansion, many modifications to General Relativity (GR), the theory 
on which modern cosmology is based, have been proposed recently. One such 
modification consists of a class of higher-order gravity models that attempt 
to address the shortcomings of GR in the infrared (IR) and ultraviolet (UV) 
ranges [1, 2, 3, 4, 5, 6, 7, 8]. These models are generally obtained by including 
higher-order curvature invariants in the Einstein-Hilbert action, by making 
the action nonlinear in the Ricci curvature i?, or contain terms involving com¬ 
binations of derivatives of R, in which case the models are known as f{R) 
theories of gravity. 

First proposed by Buchdal [9], f{R) theories gained more popularity after 
further developments by Starobinsky [10] and later following the realization of 
the discrepancy between theory and observation [11, 12, 13, 5, 14, 15]. 

The role of shear in general relativistic [16, 17, 18, 19, 20, 21, 22, 23, 24, 
25, 26] and f{R) cosmologies [27, 28, 29, 30, 31, 32, 33, 34] has been the 
subject of intense study for some time now, with the studies focusing mostly 
on the special nature of shear-free cases. In particular, it was shown in [16] 
that in the orthogonally spatially homogeneous models with vanishing shear, 
the anisotropic stresses are related to the anisotropic curvature of the spatial 
hypersurface through the electric part of the Weyl tensor. It was also shown 
that within the class of orthogonal models, small perturbations of shear are 
damped, and that the electric part of the Weyl tensor and the anisotropic 
stress tensor decay with the expansion. 

The main focus of this work is the analysis of anisotropic but homogeneous, 
shear-free models whose underlying theory of gravitational interaction is f{R)- 
gravity. 

The rest of this paper is organised as follows: in Sec. 2 a covariant descrip¬ 
tion of f{R) field equations is presented. In Sec. 3 we specialise to orthogonal 
cosmological models with anisotropic matter sources and analyse the proper¬ 
ties of such models in the case of shear-free imperfect fluids in Sec. 4. In Sec. 
5, the analysis is taken further by considering subclasses of locally rotation- 
ally symmetric spacetimes with barotropic equations of state and a qualitative 
analysis of such models has been made. Finally in Sec. 6 we discuss the results 
and give conclusions. 

Natural units [h = c = ks = 8 ttG = 1) will be used throughout this 
paper, and Latin indices run from 0 to 3. The symbols V, V and the overdot ■ 
represent the usual covariant derivative, the spatial covariant derivative, and 
differentiation with respect to cosmic time. We use the (— I- -I—b) spacetime 
signature and the Riemann tensor is defined by 


rya _ T^a r^a , T~<e r^a r^f r~ia 

^bcd ^ bd,c ^ bc,d ^ bd^ ce ^ bc^ df •> 
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where the are the Christoffel symbols (i.e., symmetric in the lower indices), 
defined by 

^bd ~ 2 ^ {9be,d F Qed.b 9bd,e] • 

The Ricci tensor is obtained by contracting the first and the third indices of 
the Riemann tensor: 

Fab — 9 Fcadb ■ 

The completely anti-symmetric pseudotensor is defined such that 

t/0123 = , 

where g = det(go&) is the determinant of the metric^ah- 

Unless otherwise stated, primes etc are shorthands for derivatives with 
respect to the Ricci scalar 

R = R\ 

and / is used as a shorthand for f{R). Moreover the following standard nota¬ 
tions are used: 

(ab) : symmetrization over the indices a and b, 

[ah\ : anti-symmetrization over the indices a and b, 

(ab) : orthogonal, symmetric, trace-free projection over the indices a and b. 


2 Covariant Description of the Field Equations 

In the standard /(i?)-gravity formulation, one starts with the modified Einstein- 
Hilbert action 

^ = i y d^Xy/^ [f{R) + 2C^] , (1) 

where Cm stands for the matter field contribution to the Lagrangian, and uses 
the variational principle of least action with respect to the metric gab to obtain 
the generalised Einstein Field Equations (EFEs) 

Gab=f^b+T^b = Tab. (2) 

Here we have defined 

-| 

= - [i(/ - RHgab + VbVaf - gabVcV^f] (3) 

as the effective matter and curvature energy-momentum tensors (EMTs), re¬ 
spectively. The EMT of standard matter is given by 

T™ = ^ = RMb + Pmhab + qTub + qTUa + , (4) 
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where /Um, Pm, 9™ and are the associated energy density, isotropic pressure, 
heat flux and anisotropic pressure, respectively, and u°‘ = is the normalized 
4-velocity of fundamental observers comoving with the fluid. We use this vector 
to define the covariant time derivative for any tensor along an observer’s 
worldlines: 


Qa..h 

^c..d 




( 5 ) 


On the other hand, we use the projection tensor hab = gab + UaUb to define 
the fully orthogonally projected covariant derivative for any tensor 


= hW..:h';,hlh:VrSH 


( 6 ) 


with total projection on all the free indices. We extract the orthogonally pro¬ 
jected symmetric trace-free part of vectors and rank-2 tensors using 





1 tab 




^cd 


( 7 ) 


and the volume element for the restspaces orthogonal to is given by [35] 

^abc — ^ Vdabc — \/\^^[a^b^c^ d]^ ^abc — ^[abc], ^abc'd’ — 0, (8) 

where rjabcd is the 4-dimensional volume element satisfying the conditions 


gabcd — Vlabcd] — ‘^^ab[c'^d] ‘^'^[a^b]cd' ( 9 ) 

The covariant spatial divergence and curl of vectors and rank-2 tensors are 
given as [36] 

divl/ = V“t4 , (div5)a = V’^Sab , (10) 

CUrlVa = Sabc^^V^" , CUrlSab = £cd(aS‘'Sb)'^ ■ (H) 

The 4-velocity vector field can be split into its irreducible parts as follows 

^a'^b — -^a'^b T "^^ab^ T ^ab “t“ ^abc^ ; 

where Aa = iia, O = WaU^-, cfab = V(aU&) and a;“ = e'^^’^'^bUc ■ 

We can also split the Weyl conformal curvature tensor [35, 37] 

C^Kd = R'^\d - + I(13) 

into its “gravito-electric” (GE) and “gravito-magnetic” (GM) parts, respec¬ 
tively, as 

Eab = CagbhU^U^, Hab = CghbdU^u'^ ■ (14) 


The GE and GM components represent the free gravitational field [35] and 
they describe gravitational action at a distance - tidal forces and gravitational 
waves. They influence the motion of matter and radiation through the geodesic 
deviation for timelike and null-vector fields, respectively. 
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The total energy density, isotropic and anisotropic pressures and heat flux 
of the f{R) universe are given, respectively, by [38] 


_ Rm 

A* = ’ 


_ Pm 

p = ~jr+PR ’ 


'^ab — 


T 


^ab ") 


a 


r,R 


(15) 


where the thermodynamic quantities for the curvature fluid component are 
defined as 




^{Rf-f)-0f"R + rv^R 
l{f-Rn + f"R + f'"R^ 


+ I (0/"i? - rV^R - fV^RVaR) , 


fR^aR + f'VaR - ^f'OVaR , 
/"V(aV6)i? + - aabRf 


(16) 

(17) 

(18) 
(19) 


In the 1-1-3 covariant decomposition [39, 40], a fundamental observer slices 
spacetime into time and space. The Bianchi and Ricci identities 


^ [aRbc]d — 0 ) (^b ^ b^a)Rc — Rabc Rd 


( 20 ) 


applied on the total fluid 4-velocity result in evolution equations - which 
propagate consistent initial data on some initial ft = to) hypersurface Sq 
uniquely along timelike congruences - and constraint equations - which restrict 
the initial data to be specified [41]. In f{R) gravity, the evolution equations 
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are given by [38] 

Am = + Pm )6 ~ ^ Qa ~ ~ ^b'^a,m ) 

PR = -iPR+ Pr)0 + , 

0 = - i©" - + 3p) + VaA“ - AaA"^ - fTaba^' + 2a;,w“ 


( 21 ) 
( 22 ) 
(23) 

C = -|0C - (Pm + Pm)A„ - VaPm - 
ia = -lOQa + ^aPR - 


r- 


~ (PR +PR)Aa — A TT^f, — £abc^ Qr , 

UJa = -§0Wa - iffabcV^’A^ + , 

^ab = —^OcJab — Eab + ^TTab + V(a^ 6 ) + ^( 0 ^ 6 ) ~ <^{aO'b)c ~ ^{a^b) ) 
Ab + ^TTab = ecd<aV^A) “ (i;ab + ^TTab) - ^ {p + p) Uab “ 

+ 3^ ('®b)c ~ ~ A(^aQb) + Scd{a ‘2A‘^+ a;°(A) + \^b)) 

Hah = —OHab — £cd{o^^E^ + \£ cd{a^‘^''^b) 


+ 3A^^b)c + ^^{aQb) — 


cd{a 




(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


Rab+<y{ab) +9c^ab-'^{aAb) - Af^aAb) -TTab- \^P - J hab = 0 , (35) 

where Rab is the Ricci tensor on 3-D spatial hypersurfaces with i? = 2/i — 
|0^ -I- 2(7^ as its corresponding (3-curvature) Ricci scalar. 


whereas the constraints read 

(C^)a := VVafc - |Va0 + £abc (vA^ + 2^0;“) + ?„ = 0 , 
(C^)ab ■= £cd{a^‘^^b)‘^ + ^(a^b) — Hab — ‘2A(^a^b) = 0 , 
(0^)a ^^Hab + (P- + P)<^a + £abc + abd {E‘^c + 

-f 3wb - i7r“^) = 0 , 

(0^)a := V^Eab + iV^ab - iV„/r + ^Oqa 

- ^aaQb - ^UJ^Hab - SabcW^^^H^ “ = 0 , 

(05) := V“wa - = 0 . 

The GauB-Codazzi equations are given by 


3 Orthogonal Models 


Following [16], the orthogonal models are characterised by the matter EMT 
representing an anisotropic fluid without heat fluxes 

Eab — PmRaRb 4 “ Pmhab 4 “ '^ab ; 


(36) 
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the matter energy density and isotropic pressure measured by an observer 
moving with the velocity In this setting, we have an irrotational and non¬ 


accelerated flow of the vector field u°‘ and therefore uja = 0 = Aa- Thus the 
corresponding evolution and constraint equations are given by 

Rm = ~{,Rm + Pm)Q ~ ^b^a,m J (3"^) 

m = -{m+ Pr)0 + , ( 38 ) 

0 = + ( 39 ) 

€ = - VaPR - , (40) 

^ab — ^ab “t“ 2^ab ^(^a^b)c i 

Eab + ^^ab = “ O {Eab + ^ + p) dab “ {aQb) 

+ , (42) 

Hab = -OHab - + 3a<^Rb}c 

+ ^£cd{a<^b)<lR > (43) 

(C*!), := |Va0 + 9f = 0, (44) 

{C*%b := ecd(aVV,)'' -Hab = 0, (45) 

{C*^)a := V^Hab + Sabc + ^bd {E'^c + = 0 , (46) 

[C*^)a := + \0q^ " 4 ^,% 

- EabcCT^'^m = 0 , (47) 


where the new equations corresponding to Eqns (26) and (34) become trivial, 
and with Eqn (24) resulting in the constraint 

(C*5)^ .= + = 0. (48) 


4 Shear-free Anisotropic Models with an Imperfect Fluid 

Erom causal relativistic thermodynamical relationships for imperfect fluids, 
the anisotropic pressure is known to evolve according to [42, 43, 44, 45] 

TTTab + TTab = -^O^ab , (49) 

where r and A are relaxation and viscosity parameters. If we consider cases 
where r is negligible and A is a positive constant, and use the fairly popular 
ansatz (valid near thermal equilibrium, such as in the very early stages of the 
Universe) for the equation of state [25, 16, 46] 

■ b ^ 


'^ab — AfJ, 


(50) 
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then Eqns (15) and (19) imply that we can rewrite (50) as 

+ /"V<aV6)i? + = aab [Rf" - A/') . (51) 

For a general case of vanishing shear tensor during the entire cosmic evolution, 
one can see from Eqn (51) that 

<b = -f"^{a^b)R - f"y^aRyb)R ■ (52) 

Moreover, the Gaufi-Codazzi equations (35) reduce to 

Rab - \Rhab =7rab = jj + f"V(aVb)R+ f"'V(aRVb}R) , (53) 

thus showing that even if the matter anisotropic stress vanishes, no constant- 
curvature geometries are guaranteed and hence no necessarily FLRW uni¬ 
verses. It is also worth noticing that, unlike in GR, if we allow the mat¬ 
ter anisotropic pressure to be nonzero despite a vanishing shear, constant- 
curvature models are allowed provided 

/"V(„Vb)R + /'"V<,i?Vb)R = 0 . (54) 

The converse also holds, i.e., it is possible, unlike in GR, to have a vanishing 
matter anisotropic pressure for a non-constant curvature geometry. 

One can see the tidal effect on the anisotropic stresses by dropping the 
shear terms of Eqn (41), obtaining the equation 

'^ab — ‘^Rab ? (55) 

which shows that, in this case as in GR [16], the anisotropic stresses are related 
to the electric part of the Weyl tensor in such a way that they balance each 
other, a necessary and sufficient condition for the shear to remain zero if 
initially vanishing. 

If the shear is nonzero, but with very small second-order contributions, 
then one can show that Eqn (41) can be approximated by 

^ab~-^0(7ab- (56) 

Rewriting Eqn (56) as 

(cr^)' Ri -|0cr2 (57) 

shows that the shear decays with expansion. One can, therefore, conclude that 
within the class of orthogonal f(R) models, small perturbations of shear are 
damped, i.e. that these models are stable if expanding, a result similar to that 
obtained in [16] for models whose underlying theory is GR. 

For shear-free orthogonal models satisfying Eqn (55), we see that Eqn (45) 
implies a purely electric Weyl tensor, i.e., Hab = 0, and hence Eqn (43) reduces 
to an identity: 

ecdia^^^E^) = ■ 


(58) 
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Moreover, it is straighforward to show using Eqns (42) and (47) that 


Eab = - ^OEab - J^iadb) > 

(59) 

V'Kb = i (V,/. - i0gf) . 

(60) 


Defining = EabE°‘^, we can rewrite Eqn (59) as 

{E^y = -lOE^ - i , (61) 

thus showing the decay of the electric part of the Weyl tensor and the anisotropic 
stress tensor with the expansion. This equation also implies decay with the 
heat flux of the curvature fluid if the bracketed terms in the r.h.s are overall 
positive. 

Let us now consider the generalized Friedman equation 

02 = 3 (/r - . (62) 

Since the total energy density /j, is not always guaranteed to be positive for 
generic f{R) models, it is not straightforward to comment on the asymptotic 
isotropization of expanding shear-free anisotropic models for the different val¬ 
ues of the spatial curvature. This is in contrast to the GR result where, for 
example, expanding shear-free models which exhibit negative spatial curvature 
asymptotically approach isotropy [16]. 


5 Anisotropic LRS Models 

Let us consider the locally rotationally symmetric (LRS) metric given by 

=—dt^ + a^{t)dr^ + b^{t) [d9^ + , ( 63 ) 

where 

{ sin(0) for R > 0 (Kantowski-Sachs), 

9 for i? = 0 (Bianchi I), 
sinh(0) for i? < 0 (Bianchi III). 

Here R = 2k/b^ for k = ±1,0. The non-vanishing kinematic quantities for 
these models are the expansion and shear, respectively given as 

e = l + ( 64 ) 

= ...a”* = ^ . (65) 

Consider the EMT of the imperfect fluid matter source to be of the form 
Tab = iJ.UaUb + phab “ 7f(ei)a(ei)b 


(66) 
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where, because of the rotational symmetry, ei = is defined as the unit 
vector along the axis of symmetry. Whereas ^ represents the total energy 
density measured by a comoving observer, the pressure measured by the same 
observer is 


P = P — TV . 


(67) 


Here the anisotropic stress tensor in the orthonormal tetrad bases 
d Id Id Id 

dt ^ a dr ’ b dO ^ hsinO d(j) 

is given by 

, / 2_ 1_ 1_\ 

TTab = diog ( 0,--7r, -TT, -TT 1 . 

This way we can write the modified EFEs as 


ab k + 


‘^ab 


62 

k + ii^ 

~1^ 


= M 


= -p + n , 


a b ab 
—^ 7 “I— 
a b ab 


( 68 ) 

(69) 

(70) 

(71) 

(72) 


whereas the conservation equations (37), (38), (40) and (48) are rewritten as 
Am — (/^m “1“ Pm 3 ^m) ^ j (^^) 


MK = - (mu +Pr- Q + 


k-mf 

f,2 


R - , 


^R _ 


C = + 


3^ya ' y/2 

VaPm(ei)“ = Va7fm(ei)“ . 


^ aR ^ aPR ^a'^R 


(74) 

(75) 

(76) 


As a result of the homogeneity assumption, tt = 7fm(t) and therefore Eqn(76) 
is trivially satisfied. 

We notice from Eqn (65) that for the case of vanishing shear, a{t) = b{t) 
and thus the modified EFES (70)-(72) reduce to 


3^ 


a 

2 “ 

a 


k 

n ^ 


a 


= -p- 


~ = ~p ■ 


Subtracting Eqn (79) from Eqn (78) yields 


(77) 

(78) 

(79) 


TT = 


k 


(80) 
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and therefore 

Eab = diag{0,-2E,E,E) , (81) 


where E = 

We adopt the barotropic EoS , Pm = (7m — l)Mm, where Pm = pm — ^ml^, 
from the continuity Eqn (73) for pm, we obtain pm = To integrate 

(77) we need to know pR. Indeed, it is a hard job to integrate (74) although we 
are working in the homogeneous case. But we can rewrite (77) in the following 
form: 



- „o „-37 v 


+ 77 — MmO 


1 

J' 





(82) 


Here p^ is the matter density at the time t = to and jm is the EoS parameter 
for the matter content. As we see, Eqn (82) is model dependent. To specify so¬ 
lutions we must choose a specific model of f{R) gravity. Otherwise, we cannot 
integrate it explicitly. Let us have a brief qualitative analysis of (82). If we are 
looking for the late-time behavior of the solutions for (82) and if we suppose 
that the space is flat k = 0, and without matter, the evolution is defined by 
the de Sitter (dS) solution, in which we put R — 6i7g, where Hg is the time 
scale of the dS universe. In this simple case, we can solve Eqn (82) to obtain: 

= (83) 

But this is not the only case we can solve (77). Suppose that we choose a 
model of f{R), so (77) reduces generally to a fourth-order ODE, which can 
be solved in terms of quadratures. For example, in the so-called Starobinsky 
model, f{R) = R + aR^, which is motivated for the inflationary universe 
scenario [10], Eqn (82) reduces to the following differential equation: 



0 -S'Vjr 

= Mmfl 


aR^ - UHR 
2 1 -k 2aR 


(84) 


where R = 6 -|- . Eqn (84 ) is a third oder ODE for a{t). So we need 

to specify initial condition(s) (ICs), as well as integrability condition(s). The 
cosmological ICs are fitted using the Hubble iJ, deceleration q, jerk j, and 
snap s parameters evaluated at the present time t = to- We can adjust the 
first derivatives of the scale factor as a(0) = og = 1 , d(0) = i7gag,d(0) = 
—iLgOg^g , d(0) = Hgjgag^ where qo is the deceleration parameter at the 
initial time (present time), jo is the jerk parameter at the instant t = to, etc 
[47]. Fortunately, these data have been measured with high precisions. 
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1.0 1i 


0 


2 


4 


Fig. 1 Numerical solution for H{t). The model of f{R) is the one proposed by Starobin- 
sky, with a = 0.02. The cosmological data are fitted with observational data for extended 
cosmological parameters. 

A series solution for a{t) in Eqn (84) has been developed using these cos- 
mographic parameters which are all evaluated at t = Iq: 


a{t) — 1 + Hq {t — to) — 1/2 Hq^qo (t — to)^ 

1 (~3 Hq^ + 54 Ho^Oi + /im + 12 a flrnHo^ 




The above solution can be used to check observational constraints. As an 
alternative, we can also solve Eqn (84) numerically. A numerical solution for 
the Hubble parameter is developed in Fig. 1 where we put ao = Hq = 1 , qq = 
—0.7. We see in Fig. 1 that H is an oscillatory function, it reaches maxima 
and minima several times. It defines an oscillatory solution but it is not in the 
form of Type IV future singularity [48, 49, 50, 51, 52]. But it can be identified 
in the late-time as the ACDM era. 

We can classify the future singularities as follows: 

— Type I: (“Big Rip”): t—a—>' 00 ,/r —>-00 and |p| —>■ 00 . 

— Type II: (“sudden”): t —>■ ts,a —>■ Og , ^ fis and \p\ —>■ 00 . 

— Type III : t —>■ tg, a —>■ Og , /r —>■ 00 and \p\ —>■ 00 

— Type IV : t —tg, a —>■ Og , /r —>■ 0 and \p\ —>■ 0 and higher derivatives of H 
diverge. Here tg , Og and ps are constants with Og ^ 0. 
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For our case, the factor given in Fig. 1, the Hubble parameter and first, second 
and third derivatives of H are plotted in Fig. 2. No higher derivatives of H 
diverges. 





-4-H . d^2H/dC2-d^3H/dt^3 

d t 


Fig. 2 Numerical solution for H , H , H- 


A phase portrait for Starobinsky model is plotted in Fig. 3. Here we solved 
the ODE with parameters 17° = 7^ = 0.3 , 7 ™ = 1. The phase portrait 

oriQ 

shows that the scale factor a(t) is a monotonic increasing function of time. It 
is always increasing, and never decreasing. 

For curiosity we are interested to know if the system has attractors or 
not. The late-time or asymptotic attractors are a class of solutions which have 
a generic form independent of the initial conditions. We examine our model 
for such types of solutions and solve the equations of motion for some initial 
conditions. The model is well established as an attractor in the following Fig. 
4. 
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6 Discussions and Conclusion 

In this work we looked at classes of shear-free anisotropic cosmological space- 
times in f{R) gravity. Focusing on orthogonal models with irrotational and 
non-accelerated fluid flows without heat fluxes, we showed that the anisotropic 
stresses are related to the electric part of the Weyl tensor in such a way that 
they balance each other. This is considered necessary and sufficient condition 
for the shear to be vanishing forever if vanishing initially. This turned out to 
be a generalization of a previous result [16] for models whose underlying the¬ 
ory is GR. We also showed that within the class of orthogonal f{R) models, 
small perturbations of shear are damped, i.e,. that these models are stable if 
expanding, and that the electric part of the Weyl tensor and the anisotropic 
stress tensor decay with the expansion as well as the heat flux of the curvature 
fluid. 

As an application, we considered a subclass of locally rotationally symmet¬ 
ric spacetimes with barotropic equations of state and studied the evolutionary 
dynamics of the Universe. In particular, we showed that the late-time be¬ 
haviour of the dS universe in f{R) gravity should satisfy Eqn (83). For the 
Starobinsky model of /(R), we provided a power-series solution for a{t) and 
we studied the behavior of the expansion parameter H (t) by numerically in¬ 
tegrating the Friedmann equation (84), where the intial conditions for Hq ,qQ 
and jo are taken from observational data. The result is the oscillatory solution 
presented in Fig. 1 and describes the late-time universe in the ACDM era. 
The first three derivatives of H have also been calculated as shown in Fig. 
2; none of these derivatives diverges. A phase-portrait anaysis for this model 
with = 0.3 ,7m = 1, given in Fig. 3, shows that the scale factor is a 
monotonically increasing function of time. Finally, we examined our model for 
late-time or asymptotic attractors, with well established solutions depicted in 
Fig. 4. 
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